clearvars;

% =========================================================================
% 设备实验数据处理，处理采集得到的超声和光声RF数据，粗重建（重建的参数应与仿真图像一致），之前项目的程序可以拿过来使用
% =========================================================================

%% define some constants to be used on the script
f0 = 2e6;       % pulse center frequency [Hz]
cycles=2;       % number of cycles in pulse
c0 = 1540;      % medium speed of sound [m/s]
rho0 = 1020;    % medium density [kg/m3]
F_number = 1.7; % F number for CPWC sequence (i.e. maximum angle)
N=1;            % number of plane waves in CPWC sequence

%% define the ultrasound probe.
%% 2.5MHZ 探头参数==================%
%% PA&US shared reconstruction parameters =========%
params.Nelements = 64;
params.fs = 40e6;   % [Hz]
params.fc = 2.5e6;  % [Hz]
params.bandwidth = 90;    % [percent]
params.pitch = 0.300e-3;  % pitch == spacing [mm]
params.width = 0.254e-3;  %  [mm]
params.radius = inf;   % 线阵
params.c = 1500;  %[m/s]
params.fnumber = [];
% xlen = params.pitch * 63 + params.width;
% x = (-xlen/2):0.05e-3:(xlen/2);
% z = 0:0.05e-3:45e-3;  % exp 60mm = 60e-3
% [X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]

% %% 10MHZ 探头参数==================%
% %% PA&US shared reconstruction parameters =========%
% params.Nelements = 64;
% params.fs = 40e6;   % [Hz]
% params.fc = 10e6;  % [Hz]
% params.bandwidth = 90;    % [percent]
% params.pitch = 0.15e-3;  % pitch == spacing [mm]
% params.width = 0.127e-3;  %  [mm]
% params.radius = inf;   % 线阵
% params.c = 1500;  %[m/s]
% params.fnumber = [];
% xlen = params.pitch * 63 + params.width;
% x = (-xlen/2):0.025e-3:(xlen/2);
% z = 0:0.025e-3:22.5e-3;  % exp 60mm = 60e-3
% [X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]

%% define the computational grid as a uff.linear_scan strcuture. 
%% We set different resolution options depending on frequency reference speed of sound.
f_max = 1.2*f0;
lambda_min = c0/f_max;

% mesh resolution, choose one
mesh_resolution='element4'; 
switch mesh_resolution
    case 'element2' % around 50 sec per wave
        dx=param.pitch/2;                            % 2 elements per pitch 
    case 'element4' % around 6min sec per wave
        dx=param.pitch/4;                            % 4 elements per pitch 
    otherwise
        error('Not a valid option');
end

% mesh size
PML_size = 20;                                          % size of the PML in grid points
Nx=round(20e-3/dx);    %% x方向的视野
Nx=Nx+mod(Nx,2);
Nz=round(40e-3/dx);    %% z方向的视野
Nz=Nz+mod(Nz,2);
grid_width = Nx*dx;
grid_depth = Nz*dx;
% domain=uff.linear_scan('x_axis', linspace(-grid_width/2,grid_width/2,Nx).', 'z_axis', linspace(0,grid_depth,Nz).');

%%%%%%%%%%%%%%%%%%%%%%%   超声重建   %%%%%%%%%%%%%%%%%%%%%%
%%%  TODO： 需要将探头收集到的超声RF数据存盘以备后续网络使用

%% Beamforming，To beamform we define a new (coarser) uff.linear_scan. We also define the processing pipeline and launch the beamformer
param.fs = 1/kgrid.dt;
% txdel = txdelay(0,100,param); % in s
xnum=200;
ynum=200;
I=zeros(200,200,N)
for n=1:3
% txdel = zeros(1,128);
txdel = delay_values(:,n);
IQ = rf2iq(element_data(:,:,n),param.fs,param.fc);
iq = IQ(:,64);
plot(abs(iq))
set(gca,'YColor','none','box','off')
xlabel('time (\mus)')
title('I/Q signal')
% legend({'in-phase','quadrature'})
axis tight

param.fnumber = [];
[xi,zi] = meshgrid(linspace(-2e-2,2e-2,xnum),linspace(0,3e-2,ynum));

bIQ = das(IQ,xi,zi,txdel,param);

I(:,:,n) = bmode(bIQ,40); % log-compressed image

end

%%%  TODO： 需要将 粗重建的超声图像存盘以备后续网络使用
figure
imagesc(xi(1,:)*1e2,zi(:,1)*1e2,sum(I,3))
colormap gray
title('PW-based echo image')
axis equal ij
set(gca,'XColor','none','box','off')
c = colorbar;
c.YTick = [0 255];
c.YTickLabel = {'-40 dB','0 dB'};
ylabel('[cm]')


%%%%%%%%%%%%%%%%%%%%%%%   光声重建   %%%%%%%%%%%%%%%%%%%%%%


%%%  TODO： 需要将探头收集到的光声声RF数据存盘以备后续网络使用
plot(element_data(:,128,1));
element_data(1:100,:)=0;
plot((element_data(:,1:14:128)/max(element_data(:))+(1:10))',...
 (0:size(element_data,1)-1)*kgrid.dt*1e6,'k')
set(gca,'XTick',1:10,'XTickLabel',int2str((1:14:128)'))
title('RF signals')
xlabel('Element number'), ylabel('time (\mus)')
xlim([0 11])
axis ij square

%% Beamforming，To beamform we define a new (coarser) uff.linear_scan. We also define the processing pipeline and launch the beamformer
param.fs = 1/kgrid.dt;
% txdel = txdelay(0,100,param); % in s
xnum=200;
ynum=200;
I=zeros(200,200,N)
for n=1:3
% txdel = zeros(1,128);
txdel = delay_values(:,n);
IQ = rf2iq(element_data(:,:,n),param.fs,param.fc);
iq = IQ(:,64);
plot(abs(iq))
set(gca,'YColor','none','box','off')
xlabel('time (\mus)')
title('I/Q signal')
% legend({'in-phase','quadrature'})
axis tight

param.fnumber = [];
[xi,zi] = meshgrid(linspace(-2e-2,2e-2,xnum),linspace(0,3e-2,ynum));

bIQ = das(IQ,xi,zi,txdel,param);

I(:,:,n) = bmode(bIQ,40); % log-compressed image

end

%%%  TODO： 需要将 粗重建的光声图像存盘以备后续网络使用

figure
imagesc(xi(1,:)*1e2,zi(:,1)*1e2,sum(I,3))
colormap gray
title('PW-based echo image')
axis equal ij
set(gca,'XColor','none','box','off')
c = colorbar;
c.YTick = [0 255];
c.YTickLabel = {'-40 dB','0 dB'};
ylabel('[cm]')




function [USPAI] = USPArecon(bin_file_path, acq_parameter_file_path, save_path, probe_param, modality, recon_method, log_compress_parameter, sound_speed)
% 这是一个提供了各种扩展可能性的US和PA图像重建的函数接口。probe_param和acq_parameter_file分别提供了探头和采集方法的参数
% Boyi RENEW 20231120 for 2.5MHz org US and PA and 10.0MHz US with EPSRI and PA
% bin_file_path            : 原始数据bin文件路径，每个bin文件第一帧为光声信号，其余帧为超声信号，此bin文件需要在前端界面选择
% acq_parameter_file_path  : 原始数据采集参数文件（.csv文件）的路径，前期用不到可忽略，无需在界面选择
% save_path                : 保存图像及图像相关参数mat文件的路径，默认保存在bin_file相同目录下，无需在界面选择
%probe_param               ：探头参数，探头参数再加上acq_parameter
% modality                 : 成像模态，0为超声，1为光声
% recon_method             : 重建算法的方法，默认为DAS，后续可能增加新的重建算法，可在界面选择
% log_compress_parameter   : 对数压缩因子，默认为14，可在界面设定
% sound_speed              : 声速，默认为1500，可在界面设定
% USPAI = USPArecon('E:\Projects\Test\Rawdata_FPGA1_Frame0.bin',1, 'E:\Projects\Test\', 0,'DAS', 16, 1500);

% disp(probe_param.Nelements);
% disp(probe_param.name);
% disp(probe_param.fc);
% disp(probe_param.pitch);
% disp(probe_param.kerf);
% disp(probe_param.width);
% disp(probe_param.bandwidth);
% disp(probe_param.radius);
% disp(probe_param.focus);
% disp(probe_param.height);




